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In an antenna array system such as that used in the Deep Space Network (DSN) 
for satellite communication, it is often necessary to account for the effects due to the 
atmosphere. Typically, the atmosphere induces amplitude and phase fluctuations 
on the transmitted downlink signal that invalidate the assumed stationarity of the 
signal model. The degree to which these perturbations affect the stationarity of 
the model depends both on parameters of the atmosphere, including wind speed 
and turbulence strength, and on parameters of the communication system, such 
as the sampling rate used. In this article, we focus on modeling the atmospheric 
phase fluctuations in a digital antenna array communications system. Based on a 
continuous-time statistical model for the atmospheric phase effects, we show how to 
obtain a related discrete-time model based on sampling the continuous-time process. 
The effects of the nonstationarity of the resulting signal model are investigated 
using the sample matrix inversion (SMI) algorithm for minimum mean-squared 
error (MMSE) equalization of the received signal. 


I. Introduction 

In a typical antenna array-based digital communications system, such as that used in the Deep Space 
Network (DSN), it is often necessary to account for effects due to the atmosphere. The atmosphere, 
which comprises the medium between the transmit antenna (typically a satellite in orbit) and the receive 
antennas (typically an array of terrestrial-based antennas), induces amplitude and phase fluctuations on 
the transmitted downlink signal [1]. These fluctuations or perturbations are often time- varying in nature, 
invalidating the stationarity of the signal model that is commonly assumed. 

If these fluctuations vary slowly enough, then commonly used adaptive algorithms for combining the 
received signals, such as the sample matrix inversion (SMI) or least mean-squares (LMS) algorithms [2], 
can be used to track the variations. On the other hand, if the atmospheric effects vary too quickly, then it 
is often difficult, if not impossible, to track the fluctuations and recover the transmitted signal. The degree 
of nonstationarity due to the effects of the atmosphere largely depends on properties of the atmosphere 
itself, including wind speed and turbulence strength, as well as properties of the digital communications 
system employed, such as the sampling rate used at the receiver. 


1 Communications Architectures and Research Section. 
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In this article, we analyze the effects of atmospheric phase variations on a digital antenna array 
communications signal. The continuous-time atmospheric phase processes are modeled by a wide sense 
stationary (WSS) statistical characterization originally introduced by Treuhaft and Lanyi [3]. From this 
continuous-time model, we derive a related discrete-time model for the atmospheric phase processes based 
on a sampled version of the phase processes. The degree of nonstationarity imparted by the atmospheric 
phase effects is analyzed experimentally through simulations using the SMI algorithm for combining of 
the array outputs. 

A. Outline 

In Section II, we review the baseband model for the antenna array digital communications signal. 
There, we show step by step the complete communication path link from the transmitter to the channel 
to the receiver. In Section III, we review the continuous-time WSS model for the atmospheric phase 
processes as derived by Treuhaft and Lanyi [3] . From this model, we show in Section IV how to obtain a 
related discrete-time WSS characterization for the sampled continuous-time atmospheric phase processes. 
In Section V, simulation results for acquisition in the presence of atmospheric phase effects are presented 
using the SMI algorithm. There, the effects due to the nonstationarity imparted by the atmosphere are 
shown. Finally, concluding remarks are made in Section VI. 


B. Notations 

All notations are as in [4] unless stated otherwise. In particular, continuous-time frequencies will be 
denoted by uppercase letters (i.e., O for radian and F for normalized frequencies), whereas discrete-time 
frequencies will be denoted by lowercase letters (i.e., u> and /). Furthermore, boldfaced lowercase letters 
(such as v) will be used to denote vectors, while boldfaced uppercase letters (such as A) will be used to 
denote matrices. One notable exception to the notation used in [4] is that continuous-time signal variables 
will be enclosed by parentheses, whereas discrete-time signal variables will be enclosed by brackets. For 
example, the notation x(t) refers to a continuous-time function of the variable t € R, whereas x[n] refers 
to a discrete-time function of the variable n G 7L. 


II. Antenna Array Digital Communications Baseband Signal Model 

Prior to introducing a statistical description of the atmospheric phase processes, it is worthwhile to 
review a fundamental model of the communications system used by the antenna array. In this section, 
we present a baseband model for the digital communications system employed by the antenna array. We 
will cover a simplified model, starting from the transmitter, progressing to the channel paths, and finally 
ending at the receiver. 

Although the goal between the transmitter and receiver is the transmission of a digital signal, the 
medium separating the two is an analog one, which requires that the digital signal be converted to an 
analog signal (such as a voltage waveform) prior to being sent through the medium. Furthermore, the 
receiver must be able to recover the desired digital signal from the received analog signals. A pictorial 
description of this scenario in the analog domain is shown in Fig. 1. 


In this setting, we have 1 transmit antenna (for example, a satellite in orbit) trying to communicate 
with N receive antennas (such as an array of terrestrial antennas). Such a system is known as a single 
input-multiple output (SIMO) system [4]. Along each path from the transmitter to the receiver, the 
transmitted signal incurs delays and distortions due to the position of the receive antennas as well as 
effects from the atmosphere [1,5]. In addition, the transmitted signal could be corrupted by noise as well 
as other interferers [5]. Prior to analyzing the effects of the analog channel, we will first review how to 
generate the transmitted analog baseband signal s(t) from a set of input digital data. 
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A. Transmitted Signal Model 

A baseband model for generating the transmitted analog signal s(t) from Fig. 1 is shown in Fig. 2. Here, 
the discrete-time digital information source d[n\ consists of complex data symbols d n . These symbols are 
typically from a special data constellation such as an M- ary phase shift keying (M-PSK) or a quadrature 
amplitude modulation (QAM) constellation [6]. 

The discrete-time signal d[n\ is converted to an analog signal using a digital-to-analog (D/A) converter. 
Here, the D/A converter generates an impulse train of Dirac delta functions [4] separated in time by an 
amount T^, which we call the baud interval or symbol interval [6]. An equivalent measure of this time is 
the baud rate or symbol rate defined by IZb = l/Tj,. The impulse train d(t) is then filtered by a pulse- 
shaping filter p(t) [6] and scaled by a positive gain factor of A to produce the output baseband transmitted 
signal s{t). 2 

Although the pulse-shaping filter p(t) can be any arbitrary shape, to prevent intersymbol interference 
(ISI) [6], the pulse typically is chosen to satisfy the Nyquist (Tf) criterion [4], given below as follows: 


TRANSMIT 

ANTENNA 



Fig. 1. Pictorial description of the analog link between the 
transmitter and receiver in an antenna array system. 


d[n] = d n 


D/A (T b ) 


d(t) = 
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Fig. 2. Baseband transmitter signal model. 


2 In a realistic implementation of the transmitter system, the impulse train generator and pulse-shaping filter operations 
are combined into one operation to generate the transmitted signal s(t) from the information source d[n]. We have opted 
to keep these operations separate here in order to make the receiver system appear more mathematically complementary 
to the transmitter system, as can be seen later on in Section II. C. 


3 


1 

Ti, 


'OO 


(1) 


/ p(t)p* (t - nT b ) dt = 6„ 


where S n is the Kronecker delta function [4] given by 


I 1, 71=0 
0, n 0, n £ 7L 


Examples of Nyquist (T b ) pulse shapes include time-limited pulses such as the non-return-to-zero (NRZ) 
pulse [6], 


Pnrz(^) — 


0 < t < T b 
otherwise 


as well as band-limited pulses such as the sine pulse [6], 


(2) 


Psinc(t) = Sine ( — 


where sine (a;) = [sin (nx) /nx]. Expressions for other common pulse shapes used, such as the Manchester, 
raised cosine, and Gaussian minimum shift keying (GMSK)-type pulses, can be found in [5,6]. 

B. Continuous-Time Channel Path Model 

Once the discrete-time information source d[n] has been converted to a continuous-time signal s(f), it 
is broadcast over the transmitter to the receiver antennas as shown in Fig. 1. Although there are many 
ways in which to model the effects of the channel, the most common way to do so is to model the channel 
as a linear time- varying distortion with additive noise [5]. This is shown in Fig. 3 for the kth channel 
path for 0 < k < N — 1. 

Here Uk(t), the output of the linear time- varying system, is given by 


Uk{t ) 


hk(r,t)s ( t - t ) dr 


(3) 


where /ife(r, t) is the impulse response of the linear time- varying system [5] [i.e., hk(t — r,t) is the out- 
put of the system in response to the input 6(t — r)]. For a terrestrial wireless setting, the impulse response 



Fig. 3. Baseband analog channel model of the Ath transmission path from Fig 1. 
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hk ( r,t ) typically consists of several multipath components along with a possible specular or line- of -sight 
(LOS) component [5]. In contrast to this, for an antenna array system such as that used in the DSN, 
the presence of multipath components will appear as very weak and the only dominant part will be the 
specular component [1]. For an array of this type, the impulse response hk (r, f) is given by 


hk{r , t) = a k (t)S(r - r k {t )) 


( 4 ) 


where ot k (t) is a complex time- varying scale factor and r k (t) is a time- varying delay parameter. 

For the DSN application under consider, the scale a k (t) accounts for both the receiver array geometry 
as well as atmospheric effects. Here, we model this scale factor as follows [1]: 


— Pk;g€ 




( 5 ) 


where p k - t9 and (p k - t9 are, respectively, the gain and phase incurred as a result of the receiver array geometry 
(both of which are assumed to be constant with time), and 9 k . a (t) denotes the phase perturbation caused 
by the atmosphere (which does vary with time). 

The delay parameter r k ( t ) typically is modeled as follows for an array system such as that used by the 
DSN: 


7~k (f ) — ^~k;g 


( 6 ) 


where r k . g is a constant delay incurred as a result of the receiver array geometry. For example, for a 
uniform linear array, the delay is given by [1] 


T k-,g — T o + 


kd 

c 


sin ip , 0 < k < N — 1 


where r 0 is a constant delay offset, d is the spacing between the array elements, c is the speed of the 
transmitted waveform in the atmospheric medium (most often the speed of light), and ’ip is the azimuthal 
angle of the transmit antenna with respect to the receive antennas [1] . 

The additive noise component rf k {t) is modeled as a circular complex additive white Gaussian noise 
(AWGN) process with mean 0 and power spectral density (psd) N 0 [5]. This is equivalent to saying that 
the real and imaginary parts of rf k {t) are uncorrelated real Gaussian processes with mean 0 and psd Nq/2. 
Mathematically, this implies that we have the following: 


E [Vk{t)] = 0 

( 7 ) 

E [Vk{t)l)k{t - t)] = 0 

(8) 

E [Vk{t)r]* k (t - r)] = N 0 S(t) 

(9) 


Furthermore, the real and imaginary parts of the processes p k (t) and rft(t) are assumed to be uncorrelated 
for all fc/f. This implies the following: 
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E[p k (to) rn (ti)] = E[r]k (to) Vi (ti)] = o ,V t 0 ,ti,fc ^ t 


(10) 


C. Received Signal Model 

Upon obtaining the analog outputs of the array channel at the receiver antennas, these continuous- 
time signals must be converted back to a discrete-time format. This typically is done using the system 
shown in Fig. 4. Note that this system is symmetric or complementary with respect to the transmitter 
system of Fig. 2. First, the received analog signal is filtered by the matched filter [6] for the original pulse 
shape pit), which is simply p* (— f). 3 The impetus for using the matched filter comes from maximum 
likelihood (ML) estimation theory [6] for estimation of the information source d[n) in the presence of 
AWGN. (It should be noted that oftentimes in an antenna-array-type environment, prior to the matched 
filtering operation, the incoming analog signal r k (t) is delay compensated [1] to account for the delays 
incurred by the channel. This compensation incorporates knowledge from predicts to estimate the path 
delays r k (t) from Eq. (6). Equivalently, this can be achieved mathematically by setting the path delays 
to zero (i.e., T k (t) = 0 for all k and t).) 



p*{-t) 

x k (t) 

a /r> (T\ 





Matched Filter Sampling Device 

Fig. 4. Baseband receiver signal model. 


j c k [n]=x k (nT) 


After matched filtering, the resulting signal x k (t) is converted to a digital signal through an analog- 
to-digital (A/D) converter [4] that simply samples x k (t) at t = nT s . Here, T s is known as the sampling 
interval [4]. Equivalently, this interval can be described by the sampling frequency or sampling rate defined 
to be F a = 1/T S . In order to recover the original information source d[n], the sampling rate F s must 
be at least equal to the baud rate 7 Zb (i.e., we need F s > Tib) [5,6]. By oversampling the output (i.e., 
when we have F s > 7 Zb), we receive additional information about the underlying transmitted signal that 
can be used to achieve a form of temporal diversity [5] . This is in addition to the inherent spatial receive 
diversity [5] that we have due to the fact that we are using multiple receive antennas. 

Assuming that the delays from the channel have been compensated and that the sampling rate F s 
equals the baud rate Tib, then from Eqs. (5), (4), and (3), and Figs. 4, 3, and 2, we have the following: 


/ OO 

i'k{r)p * (r - nT s ) dr 

-OO 

/ oo /*oo 

u k {r)p* (t - nT s ) dr + r/ k (T)p* (r - nT s ) dr 

-oo J —oo 


( 11 ) 


/ oo 

a k (T)s(r - T k {r))p* (r - nT s ) dr + v k [n\ 

-OO 


= Pk- g e^ k '' a / e 7 e fc;A T ) s ( T )p* ( T _ n T s ) dr + v k [n\ 


(12) 


3 In many digital communications textbooks, the pulse shape p(t) is time-limited to some interval t E [0,T), and the 
corresponding matched filter is implemented as p* (T — t) in order to make it causal [4]. We omit this delay here as it 
simplifies the underlying mathematics. 
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Here, Eq. (12) follows from the fact that we have assumed that the delays from the channel have been 
compensated. 

Before proceeding further with the analysis of Xk[n\, let us first consider the filtered discrete-time noise 
process Vk[n] defined in Eq. (11). Note that we have the following: 


/ OO 

E[rik{n)\p* {t - nT s ) dr = 0 

-CO 

/»oo /»oo 

E[vk[n]vk[n - m]] = / / £'[%(To)%(n)]p*(ro - nT s )p*(n - (n - m)T s ) dT 0 dri 

J — oo J — OO 


= 0 


/»oo /»oo 

- rn]] = / £[%(To)r7fc(ri)]p*(T 0 - nT s )p(r i - (n - m)T s ) dr 0 dr ± 

J — oo J —oo 

/»00 /»00 

= / / AT 0 5(to - ri)p*(r 0 - nl&)p(n - (n - m)T b ) dr 0 dri 

J — oo J — OO 

/ oo 

p*(r - nT b )p(r - (n- m)T b ) dr 

-OO 

/ OO 

p{\)p*{\-mT b )d\ 

-oo 


= 


(13) 


(14) 


(15) 


(16) 


Here, Eq. (13) follows from Eq. (7); Eq. (14) follows from Eq. (8); and Eq. (15) follows from Eq. (9) and 
the fact that T s = T b here. For Eq. (16), we used the fact that the pulse shape was assumed to satisfy 
the Nyquist (T b ) criterion of Eq. (1). Given that the continuous-time process %(t) is Gaussian, it follows 
that the discrete-time process Vk[n] is also Gaussian. From Eqs. (13), (14), and (16), we conclude that 
Vk[n] is a discrete-time AWGN process. As %(t) and %(t) satisfy Eq. (10), it can be shown that we have 


E[v k [n 0 ] ve [m] ] = E[v k [n 0 ] v\ [ni] ] = 0 , 


(17) 


and hence the discrete-time processes v bi [n] and V(\n] are uncorrelated. 

To proceed with the analysis of Xk[n], from here on in, we will assume that the pulse shape used is 
the NRZ pulse of Eq. (2). In this case, we have, from Eq. (12), 


)T S 

[n] = pk-, g e^ k ' 9 / e : ’ ek ’ a( ' T ' > s(t) dr + Vk[n] 

J nT s 


%k — Pk\g & 


(18) 


At this point, we will assume that the sampling rate F s is much larger than the rate of change of the 
atmospheric phase process. In other words, we will assume that T s is small enough that the function 
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6k;a{t) is approximately constant over the interval t £ [nT s , (n + 1)T S ). Here, we will assume that we 
have 


Ok-a{t) ~ d k . a {nT s ) , Vie [nT s , (n + 1 )T S ) 


(19) 


In Section III, we will justify this approximation. Using Eq. (19) in Eq. (18), we get the following: 


r(»+i)T. 

x k [n] = p k . g e 3 ‘ l>k ’ 9 e iek ’ a{nT ° ) / s(r) dr + v k [n] 

JnT. 


= Pfc;o e e ' 


°° /.(n+l)T» 

j<t>k-, 9p j0k ; *(nT B ) A d m I Pnrz (j ~ mT s ) dr + v k [n\ 


m =— oo 


' nT s 


= Ap k , g e^e^ nT ^ ]T d m T fe< 5 ro _„+u fc [n] 


= AT bPk , g e^e^ nT ^d n + v k [n] (20) 

Defining the discrete-time process f3 k [n] as follows, 

Pk[n] = 9 k , a {nT s ) (21) 

and using the fact that the information source d[n] is given by d[n] = d n , from Eq. (20), we have 

x k [n] = e j/3k[n] AT b p k . g e Jcl,k ' g d[n] + v k [n\ (22) 

1. Vector /Matrix Representation of the Received Array Signals. Let us define the following 
vectors and matrices: 4 


xfnl = 


xo[n\ 
xi [n] 

-X N _i[n}. 


, B[n] = diag (< 


J0o[n] p j/3 i[n] ... j0 N -i[n] 



- AT b p 0 . g e J ^°- g 


' v 0 [n\ 

A 

AThpx-ge^ 1 ’” 

r i A 

vi [n] 

a = 

- AT b p N - - 

, v[n] = 

-V N -i[n ] . 


(23) 


(24) 


Here x[n], a, and v[n] are all N x 1 column vectors whereas B[n] is an N x N matrix. From the definitions 
given in Eqs. (23) and (24), it can be seen that from Eq. (22) we have the following vector/matrix equation: 


4 Here, the notation diag (xq, xi, ■ ■ ■ ■ x ,y _ ] ) denotes an TV X N diagonal matrix whose fcth diagonal element is equal to x k . 
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x[n] = B[n]ad[n] +v[n] 


(25) 


From Eq. (25), it can be seen that x[n] consists of a linear time-varying scaled version of the desired 
scalar information source d[n\ along with an additive vector noise process v[n]. Combining Eq. (24) with 
Eqs. (13), (14), (16), and (17), it follows that v[n] is a discrete-time complex circular AWGN vector 
process [5] with mean 0 and covariance er 2 I, where we have defined a 2 as a 2 = N 0 T b . We will express 
this using the notation from [5] as 


v[n] ~ CAf(0, er 2 l) 


Note that x[n] from Eq. (25) consists of a signal component B[n]ad[n] as well as a noise component v[n]. 
One valid measure of the signal strength to the noise strength is the signal-to-noise ratio (SNR) [5,6], 
which is defined as follows: 


SNR 


x[n] 


a energy of signal component of x[n] 
energy of noise component of x[n] 


Tr 

a 

s 

PQ, 

i i 

cd 

e 

CQ, 

TT 

■e 

} 

Tr \E 

C 

3 

C 

ijr 

} 


(26) 


where Tr {•} denotes the trace operator [4] (i.e., the sum of the diagonal elements of the input matrix) and 
f denotes the transpose conjugate operator [4]. Assuming that d[n] comes from a unit energy constellation 
(so that E [|(i[n]| 2 l = 1), then from Eq. (26), we have 


SNR 


x[n] 


a i a 


Na 2 



JV-l 


N £ 


n — 0 


(27) 


The SNR as well as the combining loss (which will be defined later in Section V) will be used to measure 
the tracking performance of the SMI algorithm in the presence of atmospheric phase effects. 


III. Statistical Model for the Continuous-Time Atmospheric Phase Processes 

In the previous section, a detailed analysis of the complete communication link between the transmitter 
and receiver was given. During the course of this analysis, a set of time-varying atmospheric phase 
processes {Ok- a {t)} was introduced. In this section, we review a statistical model for this set of processes 
as based on the Treuhaft and Lanyi model given in [3] . 

Define the following N x 1 vector of atmospheric phase processes: 


m = 


do ;a(t) 

@1 ;a(t) 


L0JV— l;a(£) -I 


(28) 


Here, we model the vector process 0(f) as a wide sense stationary (WSS) random process [4]. This means 

that the mean E[6(t)] and autocorrelation E\6(t)6^(t — r)] are independent of t. As such, we can define 
the mean and autocorrelation as follows: 
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M0 = E[0(t)} 


R 0d (T)±E\o(t)0^(t-T) 


Note that even if the vector process 6{t) is WSS, the received signal vector x[n] from Eq. (25) is not 
because it contains exponentiated versions of the components of 9(t). 

Since the atmospheric phase processes are considered to be perturbations of the received signal phases, 
it is generally accepted that the atmospheric phases are zero mean, and so fJ, 0 = 0. As far as the 
second-order statistics, R^g(r), are concerned, in much of the literature, spatial correlations are modeled 
independently from temporal ones [1,3]. Due to the fact that we must jointly account for both phenomena 
here, this presents us with the dilemma of how to carry out this joint modeling. 

Before proceeding with this joint model to obtain R^g(r), we review the spatial and temporal phase 
correlation models given in [1,3]. In particular, we will focus on the models given by Treuhaft and 
Lanyi [3]. 

A. Review of Spatial and Temporal Models of the Atmospheric Phase Processes 

1. Spatial Correlation Model. The phase perturbations caused by the atmosphere are typically 
due in large part to water vapor fluctuations in the troposphere [1,3]. In these cases, the spatial refrac- 
tivity structure function [1], which we define below, can be approximated using Kolmogorov turbulence 
theory [1] . A better approximation for the spatial structure function is the one proposed by Treuhaft and 
Lanyi [3], which we describe below. 

The spatial refractivity structure function D x (r) measures the difference of the refractive index of two 
points that are a distance of r from each other. It is defined as follows [1]: 


D x (r) ± E 


lx ( r o — r ) — x ( r o)r 


(29) 


where x( r ) denotes the deviation of the refractive index from unity (i.e. , x(r) = ?y(r) — 1, where rj(r) is the 
index of refraction [3]). From Eq. (29), we have implicitly assumed that the structure function depends 
only on the distance between the two locations (ro + r) and ro, and not on the actual reference point ro- 
This stationarity assumption holds approximately for large patches of turbulence [1] . 

In classical Kolmogorov turbulence theory, D x (r) is approximated as [1] 


D x (r) = C 2 r 2/3 


where C is a parameter known as the turbulence strength [1] and is nominally 2.4 x 10 -7 m -1 / 3 . Although 
the Kolmogorov model is appropriate for small distances, it diverges as r — > oo. Physically, this does not 
make sense as the correlation between locations infinitely far away should go to zero. 

A more appropriate model for D x (r) is the one used by Treuhaft and Lanyi [3]: 


D x( r ) 


C 2 r 2 / 3 


(30) 
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where L is a parameter known as the outer scale of the turbulence [1] . Treuhaft and Lanyi chose a fairly 
large value for the outer scale (3,000 km) to be consistent with Kolmogorov’s model over a large range of 
baseline distances. More recent studies, however, conducted at the Very Large Array (VLA) in Socorro, 
New Mexico, by Carilli and Holdaway [7] indicate that a more practical value for the outer scale is 6 km. 

It can be shown that the spatial phase structure function Dg. s (r) [1], defined by 


D 9;s (r) = E 


1 0 S (r 0 + r) 


ds Ml 2 


(31) 


which measures the deviation of the phase of the atmosphere at two points separated by a distance r, can 
be given in terms of an averaged version of the refractivity structure function of Eq. (30) . In particular, 
it can be shown that we have [3] 


D, A r )= {^jc^D,Q 


(32) 


Here, A denotes the wavelength of the carrier of the transmitted signal, h denotes the height of the 
turbulence (nominally 2 km), and D s (a) is a spatial structure function given by 


D a (a) 




' a 2 + 2(x — y)a tan £ + 


x-y 
cos ^ 


- dv 



(33) 


where D x (/?) is the normalized refractivity structure function [3] given in terms of the refractivity structure 
function of Eq. (30) as 


D x m = 


i + 


/3 2 / 3 

?) 


1 

2/3 - c 2 h 2 / 3 


D x m 


(34) 


It should be noted here that £ is the zenith angle of the ray connecting the turbulence patch to the receive 
antennas as shown in Fig. 5. 

From Eq. (31), the spatial autocorrelation of the atmospheric phase Rg- S (r) is given by 

Re-,s(r) = ^[A) ;s (oo) - Dg. s {r)\ 

Using Eq. (33), Rg. s can be approximated as follows: 


Re-A r ) 


1 

2 


/ 2nChL 1 / 3 \ 2 
\ A cos £ J 


Dg- S ( r ) 


(35) 


Hence, using Eqs. (35), (32), (33), and (34), we can calculate the spatial correlations of the atmospheric 
phase processes for a fixed time. 
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A plot of Rg- S (r) is shown in Fig. 6 for the following input parameters: 

C = 2.4 x 10 -7 m -1 / 3 (turbulence strength) 
h = 2 km (height of turbulence) 

L = 6 km (outer scale) 

F c = 25 GHz (carrier frequency at Ka-band) 

Q 

A = — = 0.012 m (wavelength) 

F c 

£ = 0 (angle between turbulence patch and receive antennas measured from zenith) 

As can be seen, the phase is fairly correlated over large ranges of r. This is an expected and desired 
phenomenon, as it means that the phase perturbations seen by the receive antennas will not appear to 
be completely random. 



Receive Antennas 

Fig. 5. illustration of the position of the turbulence 
patch in relation to the receive antennas. 



Fig. 6. Spatial autocorrelation of the atmospheric 
phase R g (r) as a function of distance r. 
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2. Temporal Correlation Model. In order to derive a temporal characterization of the atmospheric 
phase processes, typically a frozen-screen assumption [1,3] is made. According to this hypothesis, the 
spatial structure of the turbulence is assumed to be frozen in time. This assumption is approximately 
true for most turbulence patches (such as clouds) and is valid over relatively long periods of time [1]. 

Assuming in addition that the turbulence patch is moving at a velocity V w (known as the wind speed) 
with respect to the receive antennas as in Fig. 5, it can be shown that we have [1,3] 


D$A T ) = D 0 . s (V w t) 


(36) 


where Dg. a (r) is the spatial phase structure function of Eq. (31) and Dg. t (r) is the temporal structure 
function defined by 


D e . t {r ) = E 


\St{T 0 + T ) 


0t(r o)| 2 


(37) 


From Eqs. (37), (36), and (31), it can be shown that we have 


Re-tir) = Rb;s(V w t) 


(38) 


where Rg-t(r) is the temporal autocorrelation of the atmospheric phase. Hence, the temporal autocorre- 
lation is simply a linearly stretched version of the spatial autocorrelation. 

Using Eq. (35), Rg- t {r) can be approximated as 


Re-A T ) 


1 

2 


/ 2irChL 1 / 3 

\ A cos £ 


2 


Dg- s {Vwr) 


(39) 


where Dg ;s (r) is given by Eqs. (32) through (34). 

A plot of Rg-t( T ) using the same input parameters as those used in Section III.A.l is shown in Fig. 7 
for various values of the wind speed V w . As expected, as V w increases, Rg-t( T ) becomes more compressed 
since the phase will appear more uncorrelated over time. In this case, this will make tracking more 
difficult since the received digital signal x[n] from Eq. (25) will change more rapidly. This difficulty in 
tracking will become more evident in Section V. 

B. Obtaining a Jointly Spatio-Temporal Characterization of the Atmospheric Phase Processes 

As stated before, although much of the literature on atmospheric phase modeling includes spatial and 
temporal statistical characterizations, there appears to be little to no work on joint characterizations. 
From this, we can infer that we can treat both effects in a separable way. In other words, we can treat 
spatial and temporal effects independently. This assumption clearly will not apply for heterogeneous 
turbulence patches, but for homogeneous ones in which the frozen-screen assumption holds, this model is 
physically plausible [1]. 

Mathematically, we model the autocorrelation function R 00( r ) separably as follows: 


R 0$( r ) = Re-A T )^e.s 


(40) 
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Fig. 7. Temporal autocorrelation of the atmospheric 
phase R g t ( r) as a function of time r. 


where Rg-t(r) is a normalized version of the temporal autocorrelation function i?e ;t (r) from Eq. (39) and 
is given by 


Re-, t( T ) 


Re-t{ T ) 

Re-A 0) 


(41) 


Here, Rg ;s is a constant matrix that depends only on the spatial autocorrelation function Rg- S (r) from 
Eq. (35). In particular, we have 


R-0;s = 


Re-A r o, o) 
Re-A r i,o) 


Re-A r o,i) 
Re-A r 1 , 1 ) 


Rg-,s{ r o,N-i) 

Re-A r i,iv-i) 


-Re-A r N- i,o) 


Re-,s( r N-i,i) 


Re-, s {i"N-i,N-i) - 


(42) 


where rk,e represents the distance between the fcth and ttYi receive antennas. As we have r^j. = r^k for 
all k. £ and ^ = 0 for all k, we can simplify Eq. (42) as follows: 


R-0;s = 


Re-A 0) 

Re-Aro, i) 


Re-A r o,i) 
Re-A 0) 


Re-A r o,N-i) 

Re-,sAi,N-i) 


-Re-A r o,jv-i) 


Re-s{ri,N-i) 


Re A 0) J 


Thus, Rg :s is a Hermitian matrix [4], as expected. For the special case in which we have a uniform linear 
array, Re/ ;s is also Toeplitz [4]. 

As can be seen from Eq. (40), R^(r) represents a model in which the spatial and temporal effects have 
been decoupled. The reason for normalizing the temporal autocorrelation as in Eq. (41) is to maintain 
consistency with the spatial autocorrelation model. Alternatively, we could have normalized the spatial 
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correlation matrix Rg :s by Rg- S ( 0) to maintain consistency with the temporal autocorrelation model. As 
Re- S ( 0) = Rg-t{ 0) from Eq. (38), both approaches yield the same expression for R 0fl( T )> as expected. 

We now proceed to derive a discrete-time statistical model for the atmospheric phase processes based 
on uniform sampling of the continuous-time vector process 0(t). There, we justify the approximation of 
the phase over a baud or symbol interval given in Eq. (19). 


IV. Discrete-Time Statistical Model for the Sampled Atmospheric Phase Processes 

A. Justifying the Sampled Phase Approximation 

Prior to introducing a discrete-time statistical description of the atmospheric phase processes based 
on sampling, we must first justify the approximation from Eq. (19) that the atmospheric phases can be 
assumed to be constant over a baud interval. In order to measure the variation of 9 k . a (t) over time, we 
can use the following quantity: 


Pfc(r 0 ,T, e) = Pr { \9 k . a (r 0 + r) - 9 k - a (T 0 )\ > e} 


(43) 


Here e is a threshold value satisfying e > 0. Note that the quantity Pfc(ro,r, e) is simply the probability 
that the magnitude of the difference of phases 9 k;a (ro + r) and 9 k;a {r o) exceeds the threshold e. 

Although it is difficult to calculate P k {r 0l r, e) from Eq. (43) in closed form, it is fairly easy to compute 
an upper bound for it using Chebyshev’s inequality [8]. Using Chebyshev’s inequality on Eq. (43), we 
have 


P k (T 0 ,T,e) < 


E 


\9 k -a(TQ + T) - 9 k - a (T 0 )\' 


De-t{r) _ D 8 . s (V w t) 


(44) 


where we have used Eqs. (37) and (36) in Eq. (44). 

Concerning Eq. (19), we are interested in studying the behavior of P k (ro,T,e) for tq = nT s and 
t € [0,T S ). From Fig. 6, it can be seen that Rg- S (r) is monotonically decreasing for r > 0. As we 
have Dg ;s (r) = 2[i?e ;s (0) — i?g ;s (r)] here by expanding Eq. (31), it follows that Dg. s (r) is monotonically 
increasing for r > 0. From Eq. (44), this implies that we have 


P k (nT s ,r,e) < 


Dg- S (V W T) Dg. s (V w T s ) 


< 


V t 6 [0, T s ) 


(45) 


For a large wind speed of V w = 50 m/s and a moderate sampling rate of F s = 1 MHz (sampling 
interval of T s = 1 fj, s), the probability that the phase changes by more than 1 degree (i.e., e = 7r/180) is 
upper bounded by 2.7535 x 10~ 4 using Eq. (45). Indeed this value is minute, and so we can say that with 
high probability the atmospheric phase is approximately constant over a sampling interval as in Eq. (19). 


B. Statistical Properties of the Sampled Phase Processes 

Recall that the discrete-time sampled phase processes /3 k [n] are defined as in Eq. (21). In this section, 
we characterize a statistical model for the discrete-time vector process /3[n\ defined below: 
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m 


Po [n] 

Pi M 

-P N -i[n\- 


Note that, from Eqs. (21) and (28), we have 


P[n] = 6{nT s ) 


Since 6(t ) is a continuous-time WSS process, it can easily be shown [4] that 0[n] is a discrete-time WSS 
process with 


fJ-ff = E[P[n ] ] = n 6 


R tf[k\ = E[0[n]fi\n - k]] = R ed (kT s ) 


(46) 


As 6{t) is assumed to be zero mean with autocorrelation as given in Eq. (40), we have 


= 0 . R pp[k] = Re-t{kT s ) R e . s 


(47) 


1. Generating Instances of the Discrete-Time Atmospheric Phase Processes. To properly 
simulate the effects of the atmospheric phase perturbations on an adaptive array algorithm as is desired 
here, we need a way to generate instances of the atmospheric phase processes. The method used to 
generate these instances is similar to the one considered in [9]. In order to accomplish this, we need to 
characterize the power spectral density (psd) of 0[n}. This quantity is simply the discrete-time Fourier 
transform of Rg^[fc] and is defined below [4]: 


OO 

Sp 0 (en = £ R 00[k)e- juk 

h — — oo 


If denotes the continuous-time Fourier transform of Rg-t(r ), defined by 


/ OO 

Re ]t (T)e- jQT dr 

-OO 


then from Eq. (47), it can be shown [4] that S^(e- Ja ’) is given by 




i £ §« 


£=—c 


tij 


u; — 2n£ 
T„ 


Re/^ 


P(e^) 


(48) 


where P(e JW ) is the aliased spectrum of i?g ;i (r) obtained due to uniform sampling 


Clip / 
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Since S$-t(jQ) is the continuous-time Fourier transform of a valid autocorrelation function it 

follows that Sg-t(j f2) is a valid continuous-time psd, meaning S$ ; t(jCl) > 0 for all O. From Eq. (48), we 
clearly have P(e JW ) > 0 for all ui, and so P(e J “) is a valid discrete-time psd, as expected. Hence, it can 
be factored as 


P{e ju> ) = \Q{e j “)\ 2 (49) 

where Q(e jaJ ) is some spectral factor [4] of P(e J “). 

Analogously, it can be shown that the Hermitian matrix R6» ;s from Eq. (42) is positive semidefinite [4], 
meaning that vtR 0;s v > 0 for all vectors v yf 0. Hence, it can be factored as follows [4]: 


Re ;s = TTt 


(50) 


Here T is some matrix spectral factor of Re ;s - 

Combining Eqs. (49) and (50) with Eq. (48), it follows that the psd S^(e J “) can be expressed as 
follows: 


Spp(e? u ) = |Q(0| 2 TTt = Q(e ju )TT^Q*(e ju ) = (' TQ{e ju )) (TQ(e J '“)) t 


(51) 


Note that Eq. (51) is a spectral factor decomposition of the psd Spp(e’ u ). 

In order to generate instances of the colored process /3[n\, we can process an uncorrelated white WSS 
process if[n\ (with = 0 and S ^{z) = I) through a spectral factor of S^(e- 7 “) as given in Eq. (51). 
If we have 


#*] = 


ipo[n] 
if i [n] 




then @[n\ can be generated as shown in Fig. 8. 

A few comments are in order here. First, as the phase perturbations are attributed to cumulative 
effects from the atmosphere, we model the distribution of the phases {(3k[n]\ to be Gaussian [1]. This 
means that the white process if[n\ leading to f3[n\ as in Fig. 8 should also be Gaussian. 

In addition to modeling the distribution of the phases, we must also consider the implementation 
of the temporal coloring filter Q(z). Here, we model Q(z ) as a stable Mth-order autoregressive (AR) 
approximation [8] to a spectral factor of P(e- JaJ ) as in Eq. (49). This means that Q(z) is implemented as 


Q(z) 


l 

W) 
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Fig. 8. System used to generate instances of the colored process pin] 
from the uncorrelated white process y/W- 

where A(z) is an Mth-order polynomial in z~ x . For the sake of stability, the zeros of A(z) are constrained 
to the region {z : |z| < 1} (i.e., inside the unit circle) [4]. The coefficients of A(z) can be found efficiently 
using the Levinson-Durbin recursion algorithm [8]. 5 

Sample realizations of the discrete-time atmospheric phase processes seen for a three-element uniform 
linear array are shown in Fig. 9 for a spacing of d = 50 m between the array elements. The atmospheric 
input parameters used were the same as those in Section III.A.l with a wind speed of V w = 10 m/s. Also, 
the order of the AR model used was M = 30. As can be seen, the phases appear to be strongly spatially 
and temporally correlated. 

The reason for such a strong temporal correlation is the high sampling rate (equivalently small sampling 
interval) used. Here, we used F s = 1 MHz, corresponding to a sampling interval of T s = 1 p,s, even though 
from Fig. 7 it can be seen that the phases typically are significantly correlated over several minutes. 

A plot of the time dependence of a sample realization of one of the phases considered above is shown 
in Fig. 10. As can be seen, the phase changes appreciably only after a long enough time (on the order of 
several seconds). This suggests that the nonstationarity caused by the atmospheric phase perturbations 
is small and should be able to be tracked sufficiently well. In the next section, this will be shown to be 
the case with respect to the SMI algorithm. 


V. Acquisition Simulation Results Using the SMI Algorithm 

Given the signals appearing at the receive antennas [equivalently the vector signal x[n] from Eq. (25)], 
we would like to obtain the desired scalar data sequence d[n\. In order to help accomplish this, the 
signals are typically combined using a linear time- varying multiple input-single output (MISO) system as 
shown in Fig. 11. Here wt [to, n] is the time-varying impulse response of the system 6 (i.e., wt [n — m,n] 
is the response to the input S n - m ) and is 1 x IV. The impetus for combining the signals in this way is 
to attempt to coherently add up the diversified versions of the desired signal d[n] as seen by the various 
receive antennas. 


5 In [9], Q(z) was modeled as a moving average (MA) approximation to a spectral factor of P(e JUJ ), meaning that Q(z) 
was itself implemented as a polynomial in z ~ 1 . This model, however, is unable to capture the non-negligible long-term 
temporal correlations that an AR model can, unless the length of the MA filter is made long enough. Making the length 
of the MA model long enough, however, is not feasible for large array simulations due to memory and computation time 
limitations. 

6 The MISO impulse response is defined using the dagger f transpose conjugate notation to make the output y\n] resemble 
the inner product of the N X 1 system w[m,n] and the shifted input x[n — m]. This is nothing more than a notational 
convention typically used in antenna array processing. 
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Fig. 9. Sample realizations of the discrete-time 
atmospheric phase processes for a three-element 
uniform linear array with a spacing of d - 50 m 
between the receive antennas: (a) first element, 
(b) second element, and (c) third element. 
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Fig. 10. Sample realization of an atmospheric phase 
process over (a) a large time scale and (b) a minute. 
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y[n]= X w ^[m,n]x[n-m\ 

m = — oo 


Fig. 11. Linear time-varying MISO system used to combine the antenna array outputs. 
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For most narrowband applications, wt [m, n\ is chosen to be a memoryless system [4] in which we have 

wt[m,n] = (^ t N, ™ = 0 =w f [n](5m 

0 , in yt 0 

In this case, we have 

y[n] = [n]x[n] (52) 

and the resulting system is memoryless in that the output is only a linear combination of the elements of 
the current input x[n] and not of past or future samples. 

A. Minimum Mean-Squared Error Combiner (Wiener Filter) 

For the purpose of acquisition, a known reference signal d r [n } often is transmitted so that the receiver 
can “phase up” to this reference. In other words, the signal combiner w[n] is chosen to make the output 
y[n] = wt[n]x[n] as close to the desired reference signal d r [n] as possible. For example, the Wiener filter 
w 0 [n] [2] minimizes the mean-squared error between the output y[n\ and reference d r [n] given by 

= E \y[n\ - d r [n\\ 2 (53) 

It can easily be shown [2] that the Wiener filter w Q [?i] can be expressed in closed form in terms of a 
matrix inverse as follows: 

w 0 [n] = K x ^ln\K xdr [n] (54) 

where we have 7 

7 l X x[fi] = E x[n]x^[?r] (55) 

U xdr [n] = E[x[n\d*[n]] (56) 

We should note that we assumed a priori that 7?. xx [?i] is invertible in arriving at Eq. (54), which in practice 
typically will be the case. If it is not, however, then the matrix inverse in Eq. (54) should be replaced 
with a pseudoinverse. 

The corresponding optimal mean-squared error £ 0 [n] is given below as follows [2]: 

= °d r M - n \ dr M^xx M (57) 

where cr^ [n] is given by a ^ [n] = E |d r [n]| . At this point, a few comments are in order. 

7 We use a calligraphic font for correlation matrices/ vectors here to distinguish the time dependence of these possibly 
nonstationary quantities from the lag dependence of WSS correlations such as the one in Eq. (46). 
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First, note that the behavior of £„[n] as a function of n can be used to gauge the joint stationarity 
of x[n] and d r [n ]. The reason for this is that if x[n] and d r [n] are jointly WSS, then £ 0 [n] is a constant. 
Heuristically, it can be argued that the more £ 0 [n] varies with n, the less jointly stationary x[?r] and d r [n } 
appear to be. We will investigate this stationarity more in detail via simulations later on in this section. 

Secondly, it should be noted that in the case in which the received signal x[n] consists of a linearly 
scaled version of d r [n] corrupted by an uncorrelated additive white noise, the Wiener filter w 0 [n] is optimal 
not only for minimizing the mean-squared error but also for maximizing the output SNR [2]. Although 
this will only be approximately true in the case in which the atmospheric phase effects are included, we 
will nonetheless still assume this to be the case when we consider the combining loss of the SMI algorithm, 
which is defined later on in this section. 

B. The Sample Matrix Inversion (SMI) Algorithm 

Although the Wiener filter is optimal for minimizing the mean-squared error of the output signal, it 
requires knowledge of the second-order statistics of the received signal x[n] (as well as the joint statistics 
of x[n] and d r [?r]). As this typically will not be known in practice, we instead estimate these statistics 
based on actual observed data. This estimation leads to the sample matrix inversion (SMI) algorithm 
that we will use here to phase up to the reference signal. 

Suppose that we have a block of K consecutive samples of x[n] from x[n — K] to x[n — 1] that were 
obtained independently of each other. To estimate the quantities 7 l K x[n] and lZ x d r [n\ from Eqs. (55) 
and (56), respectively, we use the sample mean of each of these quantities given below: 

1 K 

Kxx;K [n] = — x[n - fc]xt [n - k] (58) 

fc = i 
1 K 

'R-xdr-jcln] = — ^ x[n - k]d*[n - k] (59) 

fc = i 

If x[n — A"] , • • • , x[n — 1] and d r [n — K] , • • • , d r [n — 1] are jointly Gaussian, then the sample mean quan- 
tities of Eqs. (58) and (59) can be shown to be ML estimates of Eqs. (55) and (56), respectively [8]. 

With the sample estimates of TZ xx [n] and lZ x d r [n\ given in Eqs. (58) and (59), we can estimate the 
Wiener filter of Eq. (54) as follows: 


Here, no represents a starting time to initiate the signal combining and Kq represents the block length 
to be used initially at n = uq. 

A few comments are in order at this point. First note that, as time increases, the block length 
K = (n — no) + Kq of samples used to form the correlation estimates increases as well. Intuitively, we 
might expect that, as we collect more samples, the combining vector w[n] will approach the optimal 
Wiener filter w Q [?r]. Although this indeed is true when x[?r] and d r [n] are jointly WSS [2] (in which case 
w 0 [n] is actually a constant vector), this may not be the case when they are not. For example, if x[n] 
fluctuates rapidly in time, it may be more advantageous to fix the block size at a small value in order to 
better capture the localized variations of the signal statistics. 
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In addition to this issue, we must also exercise caution in the choice of the initial block size Kq. As 
1Z xx; k[it] from Eq. (58) is a sum of K dyadic terms (i.e., terms of the form vvt for some vector v), we 
must have K > N in order for 7 Z xx -k[it] to be invertible [4]. Alternatively, we can replace the matrix 
inverse in Eq. (60) with a pseudoinverse; however, as shown in the next section, the tracking performance 
is poor for such small values of block lengths. 

The method used to obtain the combining weight of Eq. (60) is formally the SMI algorithm. Although 
this method is computationally complex in that it requires the computation of a matrix inverse, this is 
mitigated by the fact that we can use the matrix inversion lemma [4] to obtain the new inverse given 
the previously computed one. We can also slightly reduce the complexity by removing the scaling factor 
of 1/A' present in Eqs. (58) and (59), as they cancel each other out in Eq. (60). 

C. Performance Metrics for Acquisition Evaluation 

In order to properly evaluate the acquisitioning ability of the SMI algorithm, we need a way to measure 
the performance of the SMI algorithm with respect to the mean-squared error optimal Wiener solution. 
Here, we will consider two such metrics, namely the misalignment and the combining loss, which we 
review below. 

1. Misalignment. As the Wiener filter w 0 [n] is optimal in the mean-squared sense, clearly any other 
combining weight vector w[?r] will necessarily yield at least as large a mean-squared error as w 0 [n]. One 
way to measure the performance of some combining weight with respect to the Wiener solution is to 
compute the excess mean-squared error, which is the difference between the mean-squared error of the 
combining weight under consideration and that of the Wiener filter. A normalized variant of this, known 
as the misalignment, is defined below [2]: 


r 1 A $M - go [re] _ Z[n] 

m Un] 


- 1 


(61) 


where £[n] is the mean-squared error of Eq. (53) corresponding to some combining weight vector w[n] 
and £ 0 [?r] is the mean-squared error of the Wiener filter w 0 [n] as given in Eq. (57). 

As £[n] > £ 0 [n], clearly /i[n] > 0 with equality if w[?r] is the Wiener filter w D [n]. Note that because of 
this property /i[n] is a measure of how well an adaptive algorithm can track the Wiener solution. If fj,[n ] is 
large, then the algorithm is performing poorly compared to the Wiener solution, whereas if g\n\ is small 
(i.e., close to 0), then the algorithm is performing well. For the case in which the input x[n] and reference 
signal d r [n] are jointly WSS (in which case w D [n] becomes a constant vector), the SMI and least mean 
squares (LMS) algorithms can be shown to yield /j [n] — > 0 as n — > 0 [2,10]. Furthermore, in this case, the 
misalignment of the SMI algorithm converges monotonically to 0 [10]. 

2. Combining Loss. As stated above, in the case in which x[n] consists of a linearly scaled version 
of the reference signal d r [n] along with an additive white noise component, the Wiener solution is optimal 
not only for minimizing the mean-squared error but also for maximizing the output SNR [2] . Though this 
will only approximately be true when we include the effects due to the atmospheric phases, nevertheless it 
is worthwhile to compare the output SNR of some combining weight vector w[n] with that of the Wiener 
filter w Q [n]. 

Prior to introducing such a measure, we must first properly define what is meant by output SNR in 
this setting. Consider the received signal x[n] from Eq. (25) when the reference signal d r [n\ is transmitted 
(i.e., d[n] = d r [n]). This can be expressed as 

x[n] = s[n] + v[?r] 
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where s[n] = B \n]ad r [n] is the signal component of x[n] and v[n] is the noise component. Upon combining 
the receive antenna inputs using the weight vector w[n] as in Eq. (52), we obtain the output y[n], which 
can be decomposed as 


y[n] = y s; w[n] + !/v;w[n] 


where y s;w [n] = wt[n]s[n] is the output signal component of y[n] and j/ V ;wH = wt [n]v[n] is the output 
noise component. Formally, in the spirit of Eq. (26), we define the output SNR as follows: 


SNR y;w [n] = 


E 

|j/s;w[^]| 

E 

|j/v;wM | 


With this definition of output SNR, the combining loss is simply a measure of the output SNR of some 
combining weight w[n] compared with that of the Wiener solution w G [n]. In particular, the combining 
loss C [n] is simply a measure of the difference between these two output SNRs in decibels (dB), as given 
below: 


C[n\ = 10 log 10 


/ SNR y;w [?r] \ 
V SNR y;Wo [n] ) 


101og 10 (SNR y;w [?r]) - 10 log 10 (SNR y;Wo [n]) 


= SNR y;w [n] (dB) - SNR y;Wo [n] (dB) 


(62) 


For the case in which x[n] and d r [n\ are jointly WSS and v[n] is an AWGN process, the SMI algorithm 
can be shown to have a combining loss that monotonically converges to 0 as n — » oo [10] . 

D. Simulation Results in the Presence of Atmospheric Phase Effects 

At this point, it is worthwhile to analyze the effects due to the atmospheric phase perturbations on 
the stationarity of the received signal x[n] as well as the combining performance achievable using the SMI 
algorithm for acquisition. To gauge the stationarity of x[n], we will consider the Wiener mean-squared 
error £ 0 [n] from Eq. (57), whereas to measure the combining performance of the SMI algorithm, we will 
consider both the misalignment y[n] from Eq. (61) as well as the combining loss C\n] from Eq. (62). 

For the atmospheric phases, suppose that we use the same parameters as those in Section III.A.l. 
Furthermore, suppose that the wind speed V w is set at a typical value of 10 m/s and that a 30th-order 
(i.e., M = 30) AR model was used to simulate the temporal correlation of the atmospheric phases. 

For the purpose of simulation, we opted to vary the sampling rate (equivalently the baud rate in 
this case as well) over the two values F s = 0.1 kHz and 1.0 kHz, corresponding to sampling intervals of 
T s = 10 ms and T s = 1 ms, respectively. The reason for this choice of low sampling rates is that the 
effects due to the resulting nonstationarity of the received signal caused by the atmospheric phase effects 
become more pronounced at lower rates. 

In all of the simulations considered here, we opted to use a pure tone reference coming to the antenna 
array at zenith. Here, the gain factors of both the reference and geometric array vector were set to unity 
for simplicity. In other words, we chose 
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d r [n] = 1 Vn 


a = [ 1 1 ••• 1] T 


Here, the array used was a uniform linear array with an inter-element spacing of d = 50 m. The number 
of receive antennas was chosen to be a moderate value of N = 50 here to help motivate the use of 
ultra-large arrays (ULAs) for interferometry in the DSN. Also, the gain factor A and noise variance a 2 
were chosen such that the input SNR of Eq. (26) was set to —5 dB. The reason for this choice of input 
SNR is motivated by the requirements of turbo/low-density parity-check (LDPC) codes expected to be 
implemented in the DSN array system. 


For the SMI algorithm, the starting time and initial block length were chosen to be no = 0 and Kq = 0, 
respectively. In the special case for which the block size was zero, the correlation matrices required for 
the SMI algorithm were initialized to zero. Finally, the number of ensemble realizations of x[n] for each n 
in all of the simulations was set at 8,192. 


In Fig. 12, we have plotted the observed mean-squared error £ 0 [n] for the Wiener filter as given in 
Eq. (57) for the values of F s as mentioned above. Furthermore, we have included the results obtained 
for the stationary case of no atmospheric phase effects. As can be seen, the error varied the most for the 
smallest sampling rate (largest sampling interval) and varied less for the largest sampling rate (smallest 
sampling interval). This is consistent with our intuition that the received signal x[n] should be less 
stationary the more the atmospheric phases change. The longer the sampling interval, the longer the 
atmospheric phases have to change by a significant amount from sample to sample. 


In Fig. 13, we have plotted the observed misalignment fj,[n] from Eq. (61) for each of the sampling 
rates F s from above. As before, we have included the misalignment for the stationary case of no atmo- 
spheric phase effects. From the plot, it can be seen that the misalignment in all cases reached a maximum 
value at n = 50, which corresponds to the case at which the block length K used in the SMI algorithm 
equals the number of antennas N. This is the value of the block length at which the matrix lZ xx: K[n\ 
from Eq. (58) becomes invertible. At block lengths smaller than N, is necessarily singular [4], 

whereas for block lengths greater than or equal to N, TZ xx -k [n] is most likely invertible due to the ran- 
domness of the input x[n]. The degradation seen when K = N is consistent with the theoretical results 
obtained by Reed et al. [10] for the case of no atmospheric phase perturbations. 


From Fig. 13, it also can be seen that, for n > N, the misalignment steadily decreased in all cases and 
came close to zero in all cases as well, suggesting that the SMI algorithm yielded a performance close to 
the Wiener solution. However, it can be seen that the misalignment for F s = 0.1 kHz was on average 
greater than that for F s = 1.0 kHz, suggesting that the increased nonstationarity due to the atmospheric 
phase effects made tracking the optimal Wiener solution somewhat more difficult. 


Finally, in Fig. 14, we have plotted the observed combining loss C[n] from Eq. (62). As with the 
misalignment, it can be seen that the most noticeable degradation occurred for n = N, as expected from 
the theoretical analysis of [10]. For n > N, it can be seen that the performance steadily improved and the 
loss approached zero as n increased for both sampling rates considered as well as for the stationary input 
case. This suggests that the SMI algorithm is performing like the Wiener filter in terms of output SNR 
as the block length continues to increase. As with the misalignment, the combining loss for F s = 0.1 kHz 
was on average slightly worse than for F s — 1.0 kHz, suggesting that the nonstationarity imparted by the 
atmospheric phases became more pronounced for F s = 0.1 kHz than for F s = 1.0 kHz. 
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Fig. 12. Observed mean-squared error !; o [n] for the Wiener 
filter for various sampling rates. (Results for a WSS input 
corresponding to no atmospheric phase effects included 
for comparison.) 
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Fig. 13. Observed misalignment n[n\ for various sampling 
rates. (Results for a WSS input corresponding to no 
atmospheric phase effects included for comparison.) 
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VI. Concluding Remarks 

In this article, we proposed a joint spatial/temporal statistical model for the atmospheric phase per- 
turbation processes seen at the receiver end of an antenna array system such as that used in the DSN. 
Based on this continuous-time model, we formulated a corresponding discrete-time model based on uni- 
form sampling of the continuous-time process. The impetus for considering such a description arises in a 
digital communications setting in which sampling is ubiquitously employed. 
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Fig. 14. Observed combining loss X[«] for various sampling 
rates. (Results for a WSS input corresponding to no 
atmospheric phase effects included for comparison.) 


With this discrete-time characterization in place, we analyzed the effects of the atmospheric phase 
perturbations on the acquisition performance of the SMI algorithm for combining the array outputs. 
Through simulations, it was shown that the nonstationarity imparted by the atmospheric phases was 
relatively mild for the sampling rates considered. This implies that the SMI algorithm is adequate for 
use to phase up the array outputs to a desired reference signal in these cases. 

Although the effects due to the atmosphere were shown to be minor for the moderate number of receive 
antennas used in the simulations here, it is not clear at this time how much the atmosphere will affect 
the acquisition performance when we consider a large number of antennas, such as those expected to be 
used in ultra-large arrays (ULAs). In such systems, the perturbations in phase due to the atmosphere are 
expected to be more pronounced because phases at antennas farther away from each other will be more 
weakly correlated, both in space and in time. These effects are the subject of ongoing research. 
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